function [fb_1, fb_2] = elastic_surface_force(dof_map, V, T, TE, ET, Quad_Order, FE_Type, FE_Order, bdr, fun_surface_f, varargin)
%
% Input: 
%   dof_map       ----- the dof mapping for each elements
%   (V, T, TE,ET) ----- the triangular mesh informations
%   FE_Order     ----- the order of polynomail spaces in BB Form
%   bdr               ----- the Dirichlet boundary edge list
%   g                  ----- the boundary function
%
% Output:
%   fb_1 , fb_2        ----- the force on the boundary, in two dimension
%
%  Dr. Xian-Liang Hu
%  Aug 2012
%


if(Quad_Order > 1)
    Quad_Order = 1;    % the surface force are always constants, let make it easier!
end

[qw, qx] = quad_rule_1d(6); % four point
bas_val = FE_value_1d(FE_Type, FE_Order, qx);


% some constant
cr = cr_pattern(FE_Order);
n_bdr = length(bdr);
n_dof_per_edge = (FE_Order + 1);
v_start = zeros(n_bdr,1);
v_end  = zeros(n_bdr,1);
dof_boundary = zeros(n_dof_per_edge, n_bdr);

for k = 1:n_bdr
    eg = bdr(k); tri = ET(eg,1);
    eg_local_idx = find(TE(tri,:)==eg);
    
    v_start_idx= mod(eg_local_idx, 3) + 1;
    v_end_idx = mod(v_start_idx, 3) + 1;
    v_start(k) = T(tri, v_start_idx);
    v_end(k)  =  T(tri, v_end_idx);
    
    if (eg_local_idx ~= 2)
        dof_local = cr_indices(0, FE_Order, eg_local_idx, cr);
    else
        dof_local = cr_indices(0, FE_Order, -eg_local_idx, cr);
    end
    dof_boundary(:, k) = dof_map(dof_local, tri);
    
end

% calculate the value of surface force on quadrature points
px = V(v_start, 1)*(1 - qx)' + V(v_end, 1)*qx';
py = V(v_start, 2)*(1 - qx)' + V(v_end, 2)*qx';
len = sqrt( (V(v_end,1) - V(v_start,1)).*(V(v_end,1) - V(v_start,1)) + ...
        (V(v_end,2) - V(v_start,2)).*(V(v_end,2) - V(v_start,2)) );

[g_1, g_2] = feval(fun_surface_f, px', py', varargin{:});   


n_dof = max(max(dof_map));
fb_1 = zeros(n_dof,1);
fb_2 = zeros(n_dof,1);

Jxw = len*qw';
fb_bnd_1 = (Jxw.*g_1')*bas_val;
fb_bnd_2 = (Jxw.*g_2')*bas_val;

for k = 1:n_bdr
        fb_1(dof_boundary(:,k)) = fb_1(dof_boundary(:,k)) + fb_bnd_1(k,:)';
        fb_2(dof_boundary(:,k)) = fb_2(dof_boundary(:,k)) + fb_bnd_2(k,:)';
end


end

% eval basis value on point px for order FE_Order
%
function bas_val = FE_value_1d(FE_Type, FE_Order, px)
if( strcmp(FE_Type, 'BB') )
    L1 = px;  % notice that the master element is [0,1], so that  0 <= px <= 1
    n_px = length(px);
    bas_val = zeros(n_px, FE_Order + 1);  % there are FE_Order + 1 basis
    fval = factorial(FE_Order);
    for i = 0:FE_Order
        bas_val(:, i+1) = fval/factorial(i)/factorial(FE_Order-i) * L1.^i .* (1 - L1).^(FE_Order - i);
    end
end
end